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ABSTRACT 

Recently, high resolution observations with the help of the near-infrared adaptive 
optics integral held spectrograph SINFONI (Spectrograph for INtegral Field Obser- 
vations in the Near Infrared) at the VLT proved the existence of massive and young 
nuclear star clusters in the centres of a sample of Seyfert galaxies. With the help of 
three-dimensional high resolution hydrodynamical simulations with the Pluto code, 
we follow the evolution of such clusters, especially focusing on stellar mass loss feed- 
ing gas into the ambient interstellar medium and driving turbulence. This leads to a 
vertically wide distributed clumpy or filamentary inflow of gas on large scales (tens of 
parsec), whereas a turbulent and very dense disc builds up on the parsec scale. In order 
to capture the relevant physics in the inner region, we treat this disc separately by 
viscously evolving the radial surface density distribution. This enables us to link the 
tens of parsec scale region (accessible via SINFONI observations) to the (sub-)parsec 
scale region (observable with the MIDI instrument and via water maser emission). 
Thereby, this procedure provides us with an ideal testbed for data comparison. In this 
work, we concentrate on the effects of a parametrised turbulent viscosity to generate 
angular momentum and mass transfer in the disc and additionally take star formation 
into account. Most of the input parameters are constrained by available observations 
of the nearby Seyfert 2 galaxy NGC 1068 and we discuss parameter studies for the free 
parameters. At the current age of its nuclear starburst of 250 Myr, our simulations 
yield disc sizes of the order of 0.8 to 0.9 pc, gas masses of 10 6 M and mass transfer 
rates of 0.025 Af /yr through the inner rim of the disc. This shows that our large 
scale torus model is able to approximately account for the disc size as inferred from 
interferometric observations in the mid-infrared and compares well to the extent and 
mass of a rotating disc structure as inferred from water maser observations. Several 
other observational constraints are discussed as well. 

Key words: galaxies: Seyfert - galaxies: nuclei - hydrodynamics - accretion, accretion 
discs - ISM: evolution - galaxies: structure. 



INTRODUCTION 



Often challenged, but on the whole still valid, the Unified 
Scheme of Active Galactic Nuclei ( |Antonucci|1993| |Urry fc| 
Padovani 1995 I constitutes the standard paradigm for the 



composition of Active Galactic Nuclei (AGN). Their highly 
energetic central activity is powered by accre tion onto a su- 
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permassive black hole (10 
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2 , e. g. 
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20041. Seeing the central engine or not in this scheme is 



solely determined by the orientation of a ring-like dust dis- 
tribution around the galaxy centre, the so-called molecu- 
lar torus. It constitutes not only the source of obscuration, 
but also the gas reservoir for feeding the hot accretion disc. 
Viewed face-on, most of the UV-optical continuum emission 
of the accretion disc (the so-called big blue bump) is visible in 
the spectral energy distribution (SED), whereas it is blocked 
by the dust in the case of an edge-on view. The absorbed 
energy is reemitted in the infrared wavelength regime, giv- 
ing rise to a pronounced peak in the SED of many AGN 
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( Sanders et aLp 989). This scenario enables the unification 
of two observed classes of galaxies, namely type 1 (direct 
view onto the central engine) and type 2 sources (through 
dust). It was first proposed by Antonucci & Miller (19851, 
after broad emission lines - arising from fast moving gas 
close to the central black hole (the so-called broad line re- 
gion) - have been detected in polarized light in the Seyfert 2 
galaxy NGC 1068. These photons are scattered into the line 
of sight and polarised by dust and free electrons above the 
torus opening. In unpolarised light, type 2 sources only ex- 
hibit narrow emission lines, which originate from slower gas 
motions further away from the black hole, in the so-called 
narrow line region, stretching beyond the opening of the 
torus funnel. 

Resolving the dusty torus was for the first time possible 
with the advent of the mid-infrared interferometer (MIDI) 
at the VLTI (very large telescope interferometer) by combin- 
ing the light of two of the eight meter class VLT telescopes. 
This revealed geometrically thick dust structures on parsec 



scale in a large sample of nearby AGN ( Tristram et al. 2009 1 , 
but with a large variety of torus properties. Exhaustive ob- 
servations of the nearby galaxy NGC 1068 ( jJaffe et al.|2004 



Poncelet et al. 20061 Raban et al. 2009 \ and the Circinus 



galaxy ( Tristram et al.|2007 l both found the dust to be dis- 
tributed in a complex structure. The brightness distribution 
could be disentangled into two components: a geometrically 
thin disc-like structure on sub-parsec to parsec scale and a 
fluffy or filamentary torus-like distribution surrounding it. 
For the case of the Circinus galaxy, the interferometric sig- 
nal even showed fmestructure. This was interpreted to result 
from a clumpy large scale component, which puts severe con- 
straints on theoretical models. Not being capable of direct 
imaging, radiative transfer models are needed for a careful 
data interpretation and modelling. 

From the theoretical perspective, basically three ques- 
tions are of current interest: 

O What is the structure of tori? 

© What are the basic shaping processes and how do they 
sustain the scale height of the torus against gravity? 

© How are tori built up and how do they evolve and feed gas 
on to the central engine? 

Modelling so far concentrated mainly on several aspects 
of the first two questions. A large number of continuous torus 
models have been used in radiative transfer calculations in 
order to compare simulated SEDs with high resolution ob- 



serrations in the infrared (e.g. 


Pier & Krolik||1992 


1993 


Granato & Dancse 1994) |van Bemmel & Dullemonc 


2003 


Schartmann et al. 2005 


Fritz et al. 2006|l, constraining phys- 



ical parameters of tori in specific galaxies or whole samples 
of AGN tori, like the shape, size, luminosity, etc. The most 
up to date simulations concentrate on dumpiness of the ab- 
sorbing dust (e. g. [Nenkova et al.||2002| [Honig et al.||2006 



a suggestion by|Krolik fc Begelman (1988), which helps to 



Nenkova et al.| |2008a b Scha rtmann et al.||2008 l, following 



prolong the timescale for the destruction of dust impacted 
by hot surrounding gas. Simultaneous agreement with high- 
resolution spectral (e. g. NACO or Spitzer) as well as in- 
terferometric data (MIDI) give us a good idea of possible 
structural properties of these tori. 

Concerning the second question of the structuring agent 
and the stability of the vertical scale height against gravity, 



several modelling attempts have been put forward: The torus 
was claimed to be made up of clumps with supersonic ran- 
dom velocities, which are sustained by transferring orbital 
shear energy with the help of sufficiently elastic collisions 
between the clumps (e. g. Krolik & Begelman 1988 Beck 



ert & Duschl 2004). Starburst-driven AGN tori were inves- 



tigated with the help of hydrodynamical simulations of the 
interstellar medium (ISM), concentrating on two effects: (i) 
a large rate of core-collapse supernova explosions (following 
in situ star formation in the densest region) in a rotation- 
ally supported thin disc ( Wada & Norman 2002 Wada et al 
|2009[ ) and (ii) feedback due to strong winds of young and 



massive stars ([Nayakshin fc Cuadra|[2007[ ). Alternative sce- 
narios replace the torus itself by a hydromagnetic disc wind, 
in which dusty clouds are embedded and produce the nec- 
essary obscuration differences between type 1 and type 2 



objects (Konigl & Kartje 1994 Elitzur & Shlosman 2006 \ 



or by nuclear gas discs, which are warped due to an axisym- 
metric gravitational potential generated by discs of young 
stars (e.g. |Nayakshin|2005[ ). |Pier fc Krolik| ( |1992[ ) claim the 
importance of infrared radiation pressure in the sub-parsec 
to parsec scale part of the torus or nuclear disc component. 
This scenario was elaborated in Krolik (20071. With the 
help of analytical calculations, they find good comparison 
of gas column densities to observations for reasonable in- 
trinsic AGN luminosities. For a more detailed discussion of 



available models see also Krolik (20071 



Despite these great successes in pinning down several 
aspects of physical processes in the vicinitjfjof black holes, 
a conclusive global physical model for the build-up and evo- 
lution of gas and dust structures in the nuclear region of 
active galaxies, as well as the origin of the dust content is 
still missing. The main reasons are the complexity of the 
physical processes, which are thought to happen in these re- 
gions and their scale lengths, which make them invisible for 
direct observations in most wavebands. However, a scenario 
of co-evolution of nuclear starbursts and tori, based on the 



models of Vollmer et al. (20041, 


Beckert & Duschl (2004) 


and observations of Davies et al. 


(20071 was suggested by 



three phases: (i) massive gas infall leads to a turbulent, stel- 
lar wind-driven Q w 1 disc, where Q is the Toomre stability 
parameter, (ii) subsequent supernova type II explosions re- 
move the intercloud medium and a geometrically thick col- 
lisional clumpy disc remains. In the final phase (iii) with a 
low mass accretion rate, the torus gets geometrically thin 
and transparent. Hydrodynamical simulations of the effects 
of stellar feedback from a young nuclear star cluster as ubiq- 
uitously found in Seyfert galaxies ( jDavies et al.||2007 l was 
investigated in an exploratory study by |Schartmann et al.| 
(20091. Turbulent mass input in combination with super- 



nova feedback and an optically thin cooling curve yield a 
two-component structure, made-up of a geometrically thick 
filamentary torus component on tens of parsec scale and 
a geometrically thin, but turbulent disc structure. Subse- 
quent continuum dust radiative transfer modelling results 
in a good comparison with available observational data sets. 



1 In the following, we mean the parsec scale region around the 
black hole, when we talk about the vicinity or the environment 
of the black hole! 
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Given the success of these simulations, this paper is a follow- 
up, concentrating more on the long-term behaviour of this 
scenario and justifying our assumptions with the help of 
further data comparison. This is enabled by treating the 
innermost part with a simplified effective disc model and 
by concentrating on the nearby and well-studied Seyfert 2 
galaxy NGC 1068. 

After shortly reviewing the underlying global physical 
model and its numerical realisation, the resulting state of 
the gas is discussed briefly. The gas inflow of these simu- 
lations is used in a subsequent effective treatment of the 
resulting gas surface density of the inner disc structure for a 
long-term evolution study. This enables us to link the large 
and the small scale structure in these objects. The compar- 
ison of feeding parameters as well as structural parameters 
with observations and previous modelling is done to assess 
that this model is a viable option and able to explain the 
mass assembly in the disc, the disc properties as well as the 
feeding of the central supermassive black hole. Sect, [^pro- 
vides a critical discussion, before we summarise our results 
in Sect. [5] 



2 LARGE SCALE HYDRODYNAMICAL 
TORUS MODEL 

2.1 Model and numerical realisation 

In recent years, there has been growing evidence for a direct 
link between galactic nuclear activity and star formation 
on the parsec scale ( |Davies et al.|[2007[ ). Massive nuclear 
star clusters (« 10 s Mq), as ubiquitously found in Seyfert 
galaxies significantly impact their environment during many 
stages of the evolution of their stellar population. Analyti- 



the rotation velocity drops from approximately 70 km/s to 
about 30 km/s in the tens of parsec distance range, which 
corresponds well to the values derived from the SINFONI 



cal arguments based on energy considerations (Davies et al 



2007) as well as numerical simulations (Schartmann et al 



20091 reveal that during the first phase of the evolution 



of the most massive stars towards core collapse supernova, 
most of the gas will be expelled from the region of the stel- 
lar cluster. A less violent phase follows, where the mass in- 
put is dominated by low velocity winds during the Asymp- 
totic Giant Branch (AGB) phase, after approximately 50 
Myrs. At roughly this time, |Davies et al.| ( |2007[ ) observa- 
tionally find a switch-on process of the active nucleus, in- 
terpreted as the beginning of accretion onto the black hole 
of these low-velocity stellar ejecta. This is the time when 
our simulations start. The initial condition comprises of a 
very low density, equilibrium distribution of matter for the 
given gravitational potential, which is washed out soon and 
is insignificant for the further evolution. The injection of 
mass into the model space in our simulations is mimicking 
the process of slowly expanding planetary nebulae shells: a 
mass of 0.5 Mq is evenly distributed within a radius of 1 pc 
together with the equivalent of the kinetic energy of an ex- 
panding shell of 30 km/s in form of thermal energy and with 
a bulk velocity comprised of the (constant) random velocity 
and the rotational velocity (described by a power law) of 
the underlying stellar distribution. The nuclear star cluster 
in NGC 1068 has a velocity dispersion of roughly 100 km/s 
( Davies et al.|2007 |. Therefore, random motions are an im- 
portant stabilising mechanism and the bulk rotation must 
be sub-Keplerian. For the parameters of our standard model, 



observations, which scatter around roughly 45 km/s (Davies 
et al.|2007[ ). The mass input rate is chosen according to the 
formalism described in Jungwiert et al. (20011. Gas cooling 



is calculated with the help of an effective cooling curve with 
solar metallicity, prepared with the help of the Cloudy code 
( |Plewa|[T995l see also |Schartmann et al"] |2009[ Fig. 1). The 



hydrodynamical evolution is then followed within the fixed 
gravitational potential of the nuclear stellar cluster, mod- 
elled with the help of a Plummer distribution of stars and 
the Newtonian potential of a central supermassive black hole 
with an influence radius of reH = G ^F H ~ 3.4 pc. We use 
outflow boundary conditions in radial as well as polar direc- 
tion, not allowing for inflow. Periodic boundary conditions 
are applied in the azimuthal direction, where we only model 
a 90° fraction of the whole model space. 

The hydrodynamics themselves are treated with the 
help of the Pluto code (Mig none et al.|2007 1 . Being a fully 
MPQparallelized high resolution shock capturing scheme 
with a large variety of Riemann-solvers, we consider it per- 
fectly suited to treat our physical setup. For the calculations 
shown in this paper, we used the two-shock Riemann solver 
together with a linear reconstruction method and directional 
splitting on a three-dimensional spherical coordinate system. 

A more detailed description of the underlying physical 
model can be found in Schartmann et al. (2009). However, 



mind that due to large uncertainties in the delay times and 
the rates of supernova type la explosions (see discussion in 
Schartmann et al.|2009 1 , in this paper we solely concentrate 
on the effect of turbulent mass input of the nuclear star 
cluster. The basic physical parameters of our simulations 
are summarised in Table [l] They are chosen to represent 
the nearby - and therefore well studied - Seyfert 2 galaxy 
NGC 1068. 



2.2 Resulting density and temperature 
distribution 

The resulting temporal evolution of the density distribu- 
tion within a meridional plane in a three-dimensional high- 
resolution simulation (n r = 250, ng = 133, n$ = 70) with 
the parameters given in Table [I] as well as the same cut 
through the temperature distribution is shown in Fig. [I] 
In the first time steps of the simulation, the density distri- 
bution is dominated by the original, injected mass blobs. 
Due to their large number density within the domain and 
their high randomly oriented velocities, collisions are very 
frequent. In our pure hydrodynamical setup, these collisions 
are inelastic. Therefore, the blobs dissipate a large fraction 
of their kinetic energy - which is radiated away due to opti- 
cally thin line and continuum cooling processes - and many 
of them merge to form larger complexes. As discussed in 
Sect. 2.1 the stars (and therefore also the clouds at the time 



of injection) possess sub-Keplerian rotation velocities at the 
tens of parsec distance from the centre. Therefore, the larger 
gas clouds no longer possess stable orbits, after losing a cer- 
tain amount of kinetic energy (from the random component 
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Figure 1. Snapshots of the density distribution in a meridional plane after (a) 0.2 orbits (corresponding to 7 ■ 10 4 yr) and (b) 1.7 orbits 
(approximately 6 ■ 10 s yr). The panel on the right hand side (c) shows the temperature distribution in the same meridional plane after 
1.7 orbits. 



Table 1. Parameters of our standard hydro model. 



Parameter 


Value 


Reference 


M BH 


8 • 10 6 Mq 


L03 


Af* 


2.2 ■ 10 8 M Q 


D07 


Af ini 

gas 


1.0 ■ 10 2 Mq 




Rc 


25 pc 


G03 


B T 


5pc 




Kin 


0.2 pc 




Rout 


50 pc 




<T* 


lOOkm/s 


D07 


P 


0.5 




T ini 


2.0 ■ 10 6 K 




Af n 


9.1 ■ lO -10 Mq / {yr Mq ) 


.101 


Mpn 


0.5 Mq 




r 


5/3 





Mass of the black hole (Mbh), normalisation constant of the 
stellar potential (M*), initial gas mass (M™ a ), cluster core ra- 
dius (R c ), torus radius (Rt), inner radius (Ri n ), outer radius 
(Rout), stellar velocity dispersion (it*), exponent of the angular 
momentum distribution of the stars initial gas temperature 
(Ti n i), normalised mass injection rate (M n ), mass of a single in- 
jection (Mpn) and adiabatic exponent (r). The references are: 
L03 JLodato fc Bertin|2003| l, D07 (Davies et al.|2007|l, G0 3 ( |Gal-| 
|limore fc Matthews|2003| l and J01 Jungwiert et al.|(200l). For a 
detailed description of the model, we refer to |Schartmann et al.| 
(2009) . 



mainly), and tend to stream radially inwards, where they 
find their new equilibrium at their angular momentum bar- 
rier, after suffering some angular momentum redistribution 
on their way towards the centre. 

In the course of the simulation, a two-component struc- 
ture builds up: a large-scale clumpy or filamentary compo- 



nent (mainly visible in the plots of Fig. [T| and a geomet- 
rically thin, but very dense disc on parsec scale (not vis- 
ible in the plotsQ Clumps and filaments from the outer 
component fall onto the thin disc and cause some amount 
of turbulence, triggering angular momentum transport out- 
wards and accretioij^] inwards. The corresponding tempera- 
ture distribution is depicted in Fig. [ljb) and shows comple- 
mentary behaviour with respect to the density distribution: 
the densest blobs are the coldest, whereas the tenuous gas 
in-between is the hottest. This is due to the optically thin 
gas cooling, which scales quadratically with the density dis- 
tribution. Fig.[2]shows the innermost part of this simulation 
in a cut through the equatorial plane. The main disc compo- 
nent in our simulation spans a radial range between roughly 
0.5 pc and 1.1 pc. This is in approximate agreement with a 
disc found with the help of water maser observations with an 
extent between 0.65 pc to l.lpc ( jGreenhill fc Gwinn|1997 |. 
This shows that the angular momentum distribution of the 
gas flowing towards the centre seems to be reasonable. Note 
however that the outer large scale torus part is in equilib- 
rium already, but mass still assembles in the disc component 
of our simulations. Although filling most of the volume of 
the model space, the obscuration properties are mostly given 



3 To avoid confusion: The two-component structure is what we 
called the torus in Sect^ When we talk about discs in the fol- 
lowing, we always refer to this dense nuclear disc component on 
parsec scale. The accretion disc in the immediate vicinity of the 
black hole (ranging from the marginally stable orbit up to a few 
thousand Schwarzschild radii) will be called hot inner accretion 
disc from now on. 

4 Please mind that the term accretion not necessarily refers to 
accretion onto the central black hole, but it also means simply 
mass transfer through the nuclear disc in our terminology. 
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Figure 2. Cut through the innermost part of the equatorial plane 
of our 3D hydrodynamical standard model after 6- 10 5 yr, showing 
the nuclear disc component. 



by the inner dense disc component. The outer fluffy torus 
part alone can not account for the Seyfert 1/2 dichotomy. 

However, it is still a matter of debate what the actual 
driver for angular momentum redistribution in such discs 
is. Due to this missing mechanism and the lack of feedback 
due to the formation of stars in our simulations, more and 
more gas from the large-scale torus component assembles in 
the disc, which prevents us from treating the whole Seyfert 
activity cycle - which is expected to be of the order of 10 7 
to 10 8 years - in one simulation. Apart from this, these kind 
of computations are very time consuming. 

Therefore, we use our 3D hydrodynamical simulations 
to describe the large-scale structure and attach an axisym- 
metric viscous disc model at small scales. In this paper, we 
are especially interested in comparing the resulting neutral 
hydrogen column densities, disc sizes and accretion rates 
through our inner boundary to various observations. 



Being restricted by the physics we are able to take into ac- 
count in the large scale simulations, this disc continuously 
grows in mass in our simulations, with only a very small 
amount being transported inwards. Simplified effective disc 
simulations will help us in the following to test, whether a 
stationary state can be reached and how the physical pa- 
rameters of these discs compare to observations. In this ef- 
fective disc model we calculate the viscous evolution of such 
discs including star formation and taking mass input from 
our large-scale simulations into account. Angular momen- 
tum is mediated by a turbulent viscosity (e. g. generated 
by magneto-rotational instability or a self-gravitating disc) 
and star formation happens in the densest regions of the 
disc. 



Following Pringle (19811 and Lin fc Pringle (19871, the 



viscous evolution of such a disc can be described by the 
following differential equation: 



o 



(v a Z(t,R)R 3 <°g>) 



S (1) 



This equation is derived from mass and angular mo- 
mentum conservation, coupled via the radial velocity (for 
details see Pringle 1981). E(t, R) is the surface density de- 
pending on time t and radius R within the disc. The rotation 
is Keplerian yielding the following rotation frequency: 



Q = VGMbhR" 



(2) 



with the mass of the nuclear supermassive black hole Mbh 
and G being the gravitational constant. The stellar mass 
component is neglected here, as we are well inside the in- 
fluence radius of the black hole (« 3.4 pc for NGC1068). 
We assume that the viscosity is due to some turbulent pro- 
cess and describe it with the help of a turbulent a viscosity 



(Shakura fc Sunyaev||1973 1 



,H(R), 



(3) 



where ot v is a dimensionless parameter determining the 
strength of the angular momentum transport (a„ £ [0,1]), 
c s is the speed of sound and H~(R) is the scale height of the 
disc 



3 THE EFFECTIVE NUCLEAR DISC MODEL 

3.1 Overall model and numerical realisation 

The idea behind these simulations is to use the gas supplied 
by our large scale three-dimensional Pluto torus simula- 
tions, which is driven inwards due to dissipation of turbu- 
lent cloud motions in collisions. The resulting dense clumps 
cool on small time scale, acting as a stabilising mechanism 
for the clouds and preventing them from smearing out into 
a continuous medium, which happens as soon as cooling is 
turned off. Additionally, some amount of angular momen- 
tum is transported with the help of these collisions. This 
inward motion stalls at the respective radius where the an- 
gular momentum barrier is reached for the corresponding 
gas clump and leads to the formation of a nuclear gas disc. 



H(R) = 6R, 



(4) 



with the scaling parameter S. We added the extra terms 
abbreviated by 5*. They include a source of mass due to the 
input from large scale (Ei np ut(t, R), see Sect. 3.2 1 and a sink 
due to star formation (EsF(t, R), see Sect. 3.31: 



S = E- mpxl t(t,R)-E s ?(t,R) (5) 

Equation [l] is solved with the help of MATLAB'^jpdepe 
solver. A summary of the parameters of our effective disc 
standard model in addition to the ones used for the 3D 
Pluto calculations are summarised in Table [2] 



http: / / www.mathworks.com/products / matlab / 



6 M. Schartmann et al. 



Table 2. Parameters of our standard effective disc model. 



Parameter Value 



-ft hi 

flout 

6 
T 

^start 

cluster 
^end 

cluster 



0.1 pc 
100.0 pc 
0.2 
400 K 
0.05 
50Myr 
300 Myr 
5000 



Inner radius of the domain (fli n ), outer radius of the domain 
(flout), thickness of the disc (8 = disc scale height / radius of the 
disc), gas temperature (T), a viscosity parameter (»„), age of the 
nuclear star cluster at the beginning of the simulations (tS™t er ) 
and at the end (^cluster) an< ^ resolution of the simulations (n r ). 



6x1 0* 



m 4X1 4 



2x10 4 





G.O 



1 






0.25 pc / 




X 2.50 pc 




25.0 pc 





1.0 

Time [orbits] 



Figure 3. Total gas mass within the domain of our 3D PLUTO 
simulations outside a sphere of radius 0.25 pc, 2.5 pc and 25 pc. 



3.2 Mass input from our Turbulent Torus 
Simulation 

Already after the comparatively short evolution time of our 
three-dimensional hydrodynamical simulations, the outer 
part of the simulation has reached a stable equilibrium, con- 
cerning the total mass balance. This can be seen in Fig. |3j 
where the total gas mass outside a given radius (0.25 pc, 
2.5 pc and 25.0 pc) is plotted against the time. Whereas the 
smallest radius shown lies within the dense disc structure, 
we are well outside the disc further out than 2.5 pc in ra- 
dius. Therefore, we use all gas, which is advected through 
a sphere of radius 2.5 pc for our effective disc simulations. 
This amounts to gas mass transfer of 0.14 Mq/jx at a cluster 
age of 50 Myr. The temporal evolution of the accumulated 
mass of this gas is shown as the red dashed line in Fig.[5]and 
Fig. [6] The change of slope is due to the decreasing mass in- 
put rate with time, according to the analytical prescription 
of the mass loss processes of a stellar population as described 



by Jungwiert et al. (20011. The radial location of input is 
chosen such that the angular momentum of the gas com- 
ing from the large scale hydrodynamical simulations equals 



the angular momentum of a Keplerian disc at the respective 
radius: 



j 2 

J gas, 



nput 



G (Mbh + M*(R))' 



(6) 



where j gas , input is the specific angular momentum of the gas. 

This assumes that mass inflow does not interact with 
disc material before settling into the disc plane (e. g. gas 
inflow happens from all directions and the disc is geometri- 
cally thin). Fig.[4] shows the resulting histogram of the mass 
flow through a sphere of radius 2.5 pc from the Pluto sim- 
ulation. The histogram is calculated from a slightly higher 
resolved Pluto run with a resolution of n r = 200, ng = 197 
and n<£ = 103, which was restricted to a radial range be- 
tween 2.5 pc and 50 pc. The time resolution of the gas flow 
derivation is approximately 3400 yr. For the histogram, a to- 
tal of 100 snapshots have been used between an evolution 
time of 0.7 Myr and 1.0 Myr. Colour coded is the relative fre- 
quency of a given combination of radius and surface density 
growth rate within the 100 snapshots. For each radial bin 
of this histogram, the red triangle denotes the maximum of 
the respective radial bin (the most frequent combination of 
this radius and mass input rate of the snapshots). We then 
fit these triangles with a parabola and renormalise to the 
total mass input rate. This parametrised mass input rate is 
given by 



Sinput(i, R) — 10 



-12.80-0.23 R[pc} 2 



(7) 



where the radius in the disc R is given in units of parsec 
as indicated. This function, which is also shown as the red 
parabola in Fig. |4j is used in the simulations of this paper. 
Mind that a small portion of the gas at 2.5 pc possesses more 
angular momentum than it would have in a Keplerian disc 
at the same radius. 

As already mentioned before, the stars in the nuclear 
cluster rotate with sub-Keplerian velocities. If we would use 
the mass injections of the stars within the nuclear star clus- 
ter directly and feed them into a Keplerian rotating disc, 
without hydrodynamically calculating the interaction of the 
blobs, then the yellow graph would result. However, during 
the 3D PLUTO simulations, turbulent motions within the 
large scale torus component transport angular momentum 
outwards and lead to accretion of matter and therefore to a 
further concentration of the disc towards the centre, evident 
in the different shapes of the red and yellow graph. 



3.3 Treatment of star formation 

Cold nuclear discs of this kind are known to be able to 
form stars (Paumar d et al.|200 6 Nayak shin et al.|2006 Alig 



et al. 2010). In our simulations, we convert gas into stars, 



whereever the Toomre stability parameter 



Q 



ttEG 



(8) 



is smaller than unity, where is the Keplerian orbital 
frequency, c s is the local speed of sound, S is the gas surface 
density and G is the gravitational constant. The amount of 
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x 10 



R [pc] 



Figure 4. Mass input into our effective disc model, shown as 
the growth of the gas surface density at different radii, flowing in 
through a sphere with radius of 2.5 pc from our 3D Pluto tur- 
bulent torus simulation. The input radii are set according to an 
angular momentum distribution of a Keplerian disc. Shown is the 
logarithm of a histogram for 100 timesteps with a resolution of 
approximately 3400 yr between an evolution time of 0.7 Myr and 
1.0 Myr, where the 3D hydro simulation is already in an equilib- 
rium state. The red triangles denote the maxima of each radial 
bin, the red parabola is the parametrisation used as source term 
in the effective disc simulations. The yellow graph denotes the 
growth of the gas surface density, if we directly take the angular 
momentum of the mass blobs, ejected from the stellar cluster and 
do not evolve them in a hydrodynamical simulation. 



gas removed from the computational domain is chosen ac- 



cording to the Kennicutt-Schmidt star formation law (Ken- 
nicutt|1998| >: 



E SF = 2.5- 10" 



E 



M 
yr kpc 2 



(9) 



E is as usual the gas surface density within the disc. Gas 
which forms stars is simply removed from the simulation and 
will no longer participate in the viscous evolution of the disc. 

3.4 Results of the effective disc models and 
comparison with observations 

In the course of the simulation, several physical processes 
are competing with one another. Starting with a basically 
empty model space in all our simulations, the knowledge of 
the initial condition is lost on very short timescale. The mass 
inflow from evolving stars in our large scale simulation first 
builds up the rotating disc, which is then constantly replen- 
ished over time, but with a decreasing mass input rate. The 
turbulent viscosity tries to move angular momentum out- 
ward, leading to accretion of gas through the disc towards 
the inner edge and beyond. The residual of the newly intro- 
duced gas remains in the disc and at some point, when the 
Toomre stability criterion is not fulfilled anymore, gravita- 
tional collapse can occur and stars will form, contributing to 
a stellar disc, which then coexists with the gas structure. A 
comparison between the various contributions for our stan- 
dard model (for the parameters given in Tablejl] and [2} is 
shown in Fig.[5] as a function of time. The vertical dashed 
blue line marks the estimated current age of the nuclear 
star cluster in NGC 1068 of 250 Myr ( jDavies et al][20u7| |. 
The integrated gas mass introduced into the disc simulation 
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Figure 5. Comparison of mass contributions of the various com- 
ponents of our model for an a viscosity parameter of 0.05 and 
a gas temperature of 400 K. The blue dashed line denotes the 
estimated current age of the nuclear star cluster in NGC 1068. 



is given as red dashed line. Most of this mass then gets ac- 
creted through viscous processes, as shown with the yellow 
dash-dotted line. The gas residing in the disc is given as 
black solid line. It rises steeply during the first few Million 
years. By then, gas transport through the inner boundary 
and transformation of gas into stars get significant, leading 
to a turn over of the disc mass curve. In the further evolu- 
tion, the disc mass stays approximately constant, whereas 
the slopes of the total accreted gas mass as well as the total 
mass in stars flatten due to the decreasing mass inflow rate. 

The same comparison of contributions to the total in- 
tegrated mass budget, but for the case of a lower gas tem- 
perature of only 50 K is shown in Fig. [6] When decreasing 
the gas temperature, the Toomre criterion indicates instabil- 
ity already at lower gas densities. Therefore, star formation 
now dominates over the gas mass residing in the disc af- 
ter roughly 110 Myr. It contributes the largest fraction of 
the inserted gas mass. The disc mass evolution now shows a 
visible decrease at later times. The dependence of star for- 
mation on the temperature of the gas in the disc can be seen 
directly in the comparison between Fig. [5] and Pig. [6} which 
shows a difference of roughly a factor of two between the 
stellar mass in the 400 K and the 50 K simulation. In reality, 
we expect the disc to consist of a multiphase medium, with 
regions of cold, warm and even hot gas due to gas cooling 
and a variety of heating processes within the disc, so that 
the star formation rate might be intermediate between the 
examples depicted here (e. g. |Schartmann et aT] 2009). 

In the following part of this section, we study how these 
processes compete with one another, when changing the 
(generally unknown) a viscosity parameter of the disc. In 
this study, it takes the values 0.05, 0.1 and 0.2. Observation- 
ally, |King et al.| ( |2007[ ) find for thin, fully ionised accretion 
discs a typical range of values between 0.1 and 0.4. Nu- 
merical studies of magnetohydrodynamical turbulence tend 
to derive up to an order of magnitude smaller values. For 
the case of circum nuclear discs, which are only partially 
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cluster age [Myr] 

Figure 7. Final gas density distribution at a cluster age of 
Figure 6. Comparison of mass contributions of the various com- 250 Myr. Superposed as blue dashed lines are the two limiting 

ponents of our model for an a viscosity parameter of 0.05 and a curves for the range of possible surf ace density distributio ns as 

temperature of 50 K. determined from MIDI observations jKishimoto ct al. 2009). 



ionised, angular momentum transport is reduced and de- 
pends strongly on the ionisation fraction of the gas. 

Fig. [7] shows the gas surface density distribution for the 
three a parameters: a v = 0.05 as black solid line, a„ = 0.10 
as red dashed line and a v = 0.20 as green dotted line at an 
age of the nuclear star cluster of 250 Myr. The qualitative 
difference between the parameter study is as expected: the 
smaller the viscosity (scaling proportional to a v ), the more 
mass is able to assemble in the disc. This is also reflected in 
the graph of the total gas mass in the disc (Fig. [8]). 

In order to be able to compare the sizes of the discs 
with MIDI observations, the following procedure is applied: 
MIDI detects hot dust (see Sect. [I]), which can only ex- 
ist beyond the sublimation radius, which was estimated to 



roughly to a (more or less model indepen dent) surface den- 



be roughly 0.4 pc for the case of NGC 1068 (Greenhill et al 



1996 Schartmann et al.|2005[ ). Therefore, we cut our discs at 
this radius, assume a constant gas-to-dust ratio and derive 
the half width at half maximum (HWHM) of the remain- 
ing structure. This then yields a HWHM of the dust disc 
between approximately 0.84 pc and 0.86 pc for the three a 
parameters at the current age of the nuclear star cluster 
in NGC 1068. Remarkably, this is very similar to the ex- 
tent of the inner component as seen with MIDI. Here, a 
half width at half maximum of approximately 0.7 pc is mea- 



sured (Raban et al. 20091. As already mentioned before, a 



warped disc as seen in water maser observations of the same 



size was found by Greenhill & Gwinn (19971. It extents be- 



tween 0.65 pc and 1.1 pc. Water maser emission traces warm 
(« 4 00 K), high-d ensity (n flj » 10 s - 10 10 cm~ 3 ) molecular 
gas ( Elitzur||l992[ ). However, one should note that it is still 
unclear how the structure detected in hot dust emission with 
MIDI can be connected to the very dense and clumpy disc 
assumed from water maser emission. Nevertheless, the find- 
ing in different modes of observations as well as in physically 
motivated simulations suggests a common origin. 

MIDI observations of a sample of nearby AGN suggest a 
common radial structure of AGN tori with a surface bright- 
ness distribution proportional to r -2 , which corresponds 



sity distribution proportional to r~ " (Kishimoto et al. 
|2009[ ). Line segments with the two limiting exponents (-1 
and 0) are superposed as blue dashed line segments in the 
radial range of the detection of a disc-like dust structure in 
NGC 1068 by MIDI with a size of roughly 0.7 pc in radius 



I Raban et al. 2009). In this radial range, the shape of the 



surface density distribution of our discs is well consistent 
with these slopes inferred from observations by Kishimoto 
et al.| ( |2009[ ), but are rather at the upper end. 



The time evolution of the total gas mass within the disc 
(Fig. [8jl shows a steep rise at the beginning of the simulation, 
until the gas disc is build up and a maximum is reached. 
Then, it is followed by a shallow decay, due to the decreasing 
gas mass input rate. Values are of the order of 2 ■ 10 6 M@ in 



the maximum. In good comparison to that, Kumar ( 1999 1 
found a disc mass of the order of 10 6 Mq with the help of a 
clumpy disc model, which accounts for the observed maser 
emission. Furthermore, a highly inclined thin, and clumpy 
ring or nuclear disc with a molecular mass of 1.3 • 1O 6 M 



Montero-Castano et al. 2009 ) has been found in our own 



galactic centre in a distance ranging from approximately 2 
to 5pc ( |Genzel et aL]|1985| |Guesten et~aT||1987[ >. This is 
particularly interesting, bec ause the black hole mass in the 
galac tic centre of 4 • 10 6 M Q jGhez et al.|2005| |Schodel et~aL 
|2009[ ) is very similar to the one in NGC 1068 and the mass 
of its circum nuclear disc is comparable to the results of our 
infall model. 

Fig. [9] shows the total accreted mass of the parame- 
ter study. As expected, larger ov-values lead to a larger 
amount of angular momentum transfer outwards and to a 
larger amount of gas accretion through the inner boundary, 
which can also be seen in the current mass transfer rate 



through the inner boundary (Fig. 10 1 in solar masses per 



year as a function of time. After an evolution time of 250 Myr 
after the cluster formation, a typical accretion rate of ap- 
proximately 0.025 Mq/yx results. It is generally unknown, 
which fraction of this mass transported through a sphere of 
radius 0.1 pc really reaches the black hole. Assuming that 
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Figure 8. Total gas mass of the disc. 
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Figure 9. Total accreted gas mass through the inner boundary 
of the domain. 



the efficiency is somewhere in between 10% and 100%, ac- 
cretion rates onto the black hole of between approximately 
0.0025 M Q /yr and 0.025 M /yr result, which are well con- 
sistent with typically observed values for Seyfert galax ies 
of a few times 10" 3 M Q /yr to 10~ 2 M Q /yr ( |Jogee||2006| ). If 
we assume that all of the gas transferred through the inner 
boundary of our model space will reach the black hole and 
that the efficiency for mass energy conversion is 10%, this 
yields a source luminosity of 1.5 • 10 44 erg/s corresponding 
to 15% of the Eddington luminosity. For radiative transfer 
calculations through a simple two-dimensional axisymmet- 
ric continuous torus model, Schartmann et al. ( 2005[ ) needed 
a value of 20% in order to obtain a good adaptation to high 
spatial resolution spectral energy data. In radiative transfer 
calculations, the source luminosity sets the scaling of the re- 
processed radiation. Determining the source luminosity is a 
difficult task for the case of NGC 1068, being a heavily ob- 
scured Seyfert 2 galaxy. Therefore, it can only be estimated 
with the help of the accretion disc radiation reflected by 
electrons above the dust distribution and is therefore geom- 



etry dependent. Pier et al. (19941 did a careful analysis of 



the geometry dependence and came up with a best estimate 
of 



9.4 ■ 10 



3.6- 10 



/rcfl 

0.01 

i erg 
s 



D 



14.4 Mpc 



(10) 
(11) 



which corresponds to 35% of the Eddington luminos- 
ity, including a quite large uncertainty However, concern- 
ing the current mass accretion rate, our smooth disc de- 
scription with an effective a parameter is an oversimplifi- 
cation. In reality, we expect the matter transport from the 
circum nuclear disc towards the inner hot accretion disc to 
be clumpy, as e. g. observed in the galactic centre environ- 



ment (Montero-Castano et al. 20091 



6 Equ. |lO| has been rescaled to today's distance estimate, com- 
pared to the original publication. 
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Figure 10. Total accretion rate of gas through the inner bound- 
ary of the domain. 



The total mass of formed stars for the a parameter 
study is shown in Fig. |11| The star formation rate increases 
steeply at the beginning of the simulation, when the mass 
input rate is still high. Larger values of a lead to a larger 



mass accretion rate (see Fig. 10 1 and a smaller gas surface 
density (see Fig. [T| and in consequence to a smaller total 
mass in stars. A total mass of up to 4 ■ 10 6 Mq is reached 
at the current age of the nuclear starcluster. In compari- 
son to this, [Kumar (1999) estimate the stellar mass within 
1 pc around the nuclear black hole in NGC 1068 to be of the 
order of 10 7 Mq on basis of theoretical considerations and 
from observations by |Thatte et al.| ( |1997[ ) . Such a high rate 
of star formation will also lead to subsequent energy input 
into the disc by supernova type II explosions after a few mil- 
lion years. This additional feedback process will stir further 
turbulence in the disc and might be an important driver 
to increase its scale height. However, a detailed analysis of 
these effects can only be done in three-dimensional, multi- 
phase hydrodynamical simulations including the effect of gas 
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Figure 11. Total mass of formed stars. 



cooling. However, using such a scenario, |Wada et al.] ( [2009] ) 
find that only on scales beyond a few parsec, a significant 
scale height can be achieved. 

A further comparison with observations can be done 
by calculating the neutral hydrogen column density on the 
line of sight. To do this, we use the gas surface density and 
distribute the gas homogeneously in z-direction in a disc 
structure according to the scale height distribution assumed 
previously (equation |4j. Many simplifications have to be as- 
sumed, concerning the geometrical distribution, the gas mix- 
ture, the ionisation state, etc. Therefore, only a rough esti- 
mate can be given. Assuming solar gas composition and that 
the disc is completely neutral, values between 4 x 10 25 cm -2 

" 1 — 0.05 result. In a 



for a v = 0.2 and 1 x 



for a v 



compilation of observed Seyfert galaxies in Shi et al. (20061, 
the values spread between approximately 10^ u cm -2 (which 
is basically the galactic foreground value) and 10 25 cm~ 2 . 
Seyfert 2 galaxies fill the upper part of this range start- 
ing at roughly 10 22 cm -2 . This shows that our derived val- 
ues lie in the upper range or even beyond the sample of 
Seyfert 2 galaxies. One should however note that a density 
stratification in vertical direction of the disc is expected. 



In particular, Mulchaey et al. ( 1992 1 find that the torus in 



NGC 1068 has a very high gas column density of the order of 
iVn = 10 25 cm -2 . An additional component on tens of par- 
sec scale, as revealed by SINFO NI, adds a column de nsity 
of approximately 8.0 x 10 24 cm" 2 ( |Sanchez et al.|2009| ). The 



hydrogen column density through our outer torus compo- 
nent (in our 3D hydrodynamical models) is only a negligible 
fraction of these values. 



4 DISCUSSION 

4.1 Keeping tori stable against gravity 

One major subject in theoretical AGN torus research is the 
question of stability of the vertical structure. A geometri- 
cally thick torus, which is stabilised by Keplerian rotation 
will soon collapse to a thin disc, when the thermal pressure 
is reduced by gas cooling or the turbulent pressure is dissi- 



pated, e. g. in collisions. Being a fast process happening on 
the order of only a few orbital periods, this contrasts with 
observations of geometrically thick structures. A short sum- 
mary of the models proposed to circumvent this problem is 
given in Sect, [I] 

In our three-dimensional hydrodynamical model, geo- 
metrically thick gas and dust structures naturally result 
from the fact that the emitting stars are organised in a verti- 
cally extended structure ( Davies et al.|20 07). These clusters 
are hot stellar systems, which are stabilised mainly by ran- 
dom motions of the stars, leading to a turbulent pressure and 
rotation is only of minor importance. The turbulent motions 
lead only to a weak stabilisation against gravity concerning 
the emitted g the blobs merge and dissipate their ran- 
dom kinetic energy component on a short timescale and are 
able to move radially inward, until they find their new equi- 
librium radius in the Keplerian disc. 

However, one should note again that the outer torus 
part of our model can not account for the Seyfert 1/2 di- 
chotomy, due to its low gas column density. This can only 
be done by the inner nuclear and very dense disc component. 
The crucial processes for puffing up this structure are still a 
matter of debate. Every process which is able to stir turbu- 
lence is a possible candidate. Turbulent motions within the 
disc will do both, puff-up the disc structure to form geomet- 
rically thick tori and drive accretion as it also leads to an 
effective viscosity. Some ideas are discussed in Sect. |4.2| 



4.2 Origin of the turbulent viscosity 

Evidence is growing that magnetic fields are most impor- 
tant to mediate redistribution of angular momentum in thin 
and ionised accretion discs. This basic idea was coined al- 
ready by |Shakura fc Sunyaev| ( |1973[ ) after setting up their 
parametrised thin disc model (Shakura-Sunyaev disc, a- 
disc, standard disc). The first numerical realisation of this 
Magneto-Rotational Instability was given by Balbus & Haw- 
ley ( 1991 1 showing that a weak magnetic field together with 



a shear flow is able to maintain a magnetic dynamo in ac- 
cretion discs and lead to an accretion flow. This scenario 
still forms the basis of the most up to date simulations. In 
our case, this might be relevant only for the innermost re- 
gion of the disc and a boundary layer, which can be directly 
ionised by the central source (e. g. |Blaes fc Balbus||1994[ |. 
However, there is still contradiction between theoretically 
derived a-values and those measured in observations ( |King| 
et al.|2007 l. 

Furthermore, any kind of turbulent process will lead 
to an effective viscosity and thereby drives angular momen- 
tum redistribution. Relevant for our simulations is first the 
merging of blobs and streams of gas into the nuclear disc. 
Possessing a range of various momenta coming from differ- 
ent directions, this leads to enhanced turbulence within the 
disc. Detailed hydrodynamical simulations of these processes 
are planned and will give us an estimate of the relevance for 
angular momentum transfer in the disc. 



|Rice fc Arm itagc ( 2009 ) proposed self-gravity to be the 
dominant transport mechanism of angular momentum for 
the case of discs, which are too weakly ionised to sustain 
MHD turbulence ( |Blaes fc Balbus|1994[ ). This might be the 
situation within some radial extent of cold and dense proto- 
planetary accretion discs and might be relevant for the case 
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of our nuclear discs as well. Furthermore, it has been shown 



that fragmentation is equivalent to an a ^ 0.06 (Rice et al 



2005) 



Due to the large amount of star formation in some of 
our simulations, secondary feedback processes like supernova 
type II explosions might contribute significantly during some 
phases of the evolution. Interaction with wind and ultra- 
violet radiation pressure will play a role in the part of the 
central region with direct lines of sight towards the nucleus. 



4.3 Dependence on the cluster rotation profile 

For the case of this paper, a simple power law is assumed 
for the radial dependence of the rotation velocity of the un- 
derlying nuclear stellar cluster. On the tens of parsec scale, 
this is a fair representation of the flat rotation curve ob- 



served with SINFONI (Davies et al. 20071 for NGC 1068 



Given the good comparison of our results with observations 
as discussed in Sect. |3.4| we think that our model is a pos- 
sible solution. A different rotation velocity distribution of 
the stars - as has been observed for other nearby Seyfert 
galaxies - will lead to a change of the radial dependence of 
the mass input rate, used as source term for the effective 
disc simulations. Test calculations with several analytic in- 
put rates showed only slightly different results. However, a 
detailed analysis of various dynamical initial conditions ne- 
cessitates an extensive parameter study and a detailed un- 
derstanding of the angular momentum redistribution, which 
is beyond the scope of the current paper and will be dis- 
cussed in subsequent publications. It will finally enable us 
to make statistical statements and compare our results to 
observations of samples of galaxies. 



5 CONCLUSIONS 

In this paper, we show that evolving stars from a massive 
and young nuclear star cluster, as found in nearby Seyfert 
galaxies provide enough gas to assemble a parsec-sized nu- 
clear gas disc. To this end, we combine an enhanced ver- 



sion of the Schartmann et al. ( 2009 ) models with a simpli- 



fied treatment of the innermost parsec scale region, where 
a nuclear disc builds up. As far as possible, we derive input 
parameters of our models from observations of the nearby 
and well-studied Seyfert 2 galaxy NGC 1068. This two-stage 
analysis enables us to (i) do a long term evolution study, 

(ii) link the tens of parsec scale region of galactic nuclei 
(observed with the SINFONI instrument) to the sub-parsec 
scales (probed by MIDI and in water maser emission) and 

(iii) test our model directly with a large number of obser- 
vational results. Such comparisons show that the total mass 
and surface density distribution of these nuclear discs are 
compatible with recent observations. 

Our three-dimensional hydrodynamic simulations are based 
upon the observations of young and massive nuclear star 
clusters in the centres of nearby Seyfert galaxies. In our 
global scenario, the gas ejection of their stars provide both, 
material for the obscuration within the so-called unified 
scheme of active galactic nuclei and a reservoir to fuel the 
central, active region. Gas is injected in form of blobs formed 
by emissions of single intermediate mass stars with low ex- 
pansion velocities. In contrast to this, turbulent velocities 



as taken over from the hot stellar system are high and con- 
tribute significantly to the clumping of the gas distribution 
as a whole, through interaction between gas blobs, trans- 
forming kinetic energy into heat. Cooling instability further 
confines these clouds and filaments, which then tend to move 
towards the minimum of the effective potential. After only 
a few orbital periods, this represents a clumpy and filamen- 
tary flow of gas towards the inner region, where it opens out 
into a disc structure and sets into a stationary state in the 
large scale part. As we currently cannot treat the physics of 
the inner part correctly in our 3D hydrodynamical simula- 
tions, we feed the mass inflow into our viscous disc model, 
in order to be able to estimate observable quantities. At the 
current age of the nuclear starburst in NGC 1068 of 250 yr, 
our simulations yield disc sizes of the order of 0.8 to 0.9 pc, 
which compares well to the observed maser disc extent and 
disc sizes as inferred from interferometric observations in the 
mid-infrared. The gas mass of a few times 10 6 Mq in the disc 
is in good comparison to models of discs, tuned to reproduce 
these maser observations. The resulting mass transfer rate of 
0.025 Mq /yr through the inner rim of our disc compares well 
to typical Seyfert galaxies. However, NGC 1068 seems to be 
in a heavily accreting state, which can be accomodated in 
our model only, when assuming e. g. clumpy accretion. Fur- 
thermore, rough estimates of gas column densities of our 
model are in agreement with NGC 1068 being a Compton 
thick source. On basis of these comparisons, we conclude 
that the proposed scenario seems to be a reasonable model 
and shows that nuclear star formation and subsequent AGN 
activity are intimately related. 
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